package scu.maqiang.matlabcases;

import scu.maqiang.elasticity.Elasticity3DT42;
import scu.maqiang.fes.BVPType;
import scu.maqiang.mesh.*;
import scu.maqiang.numeric.Direct;

public class Bracket {

	public static void main(String[] args) {
		String meshFile = "matlabMeshes\\Bracket.matmsh";
		Mesh3T4 mesh3D = new Mesh3T4().readMesh(meshFile, MeshFormat.matMsh);
		mesh3D.toTecplot("matlabMeshes\\Bracket.dat");
		ScalarFunc fixBd = (xy, label, param) -> Math.abs(xy[0]) < 1.0e-6?1.0:-1.0;
		ScalarFunc pressBd = (xy, label, param) -> Math.abs(xy[0] - 0.215)< 1.0e-6?1.0:-1.0;
		int fixBdLabel = 4;
		int pressBdLabel = 8;
		mesh3D.setBoundaryLabel(fixBd, null, fixBdLabel);
		mesh3D.setBoundaryLabel(pressBd, null, pressBdLabel);
		
		Elasticity3DT42 e = new Elasticity3DT42(mesh3D);
		e.elastoStaticProblem(BVPType.COMMON, new double[] { 2.0e11, 0.3 }, null,
				             new int[][] {new int[] {fixBdLabel}}, new Direct[] {Direct.All}, new double[][] {new double[] {0.0, 0.0, 0.0}}, 
				             new int[] {pressBdLabel}, new double[] {0.0, 0.0, -10000.0}, "matlabResults\\BracketDeformation.dat");
		/*
		FES3T42 fs = new FES3T42(mesh3D);
		SRMatrix A = new SRMatrix(fs.GetNdof());
		double[] RHS = new double[fs.GetNdof()];
		double[] x = new double[fs.GetNdof()];
		fs.assembleElasticStiff(new double[] {2.0e11, 0.3}, BVPType.COMMON, A);
		fs.assembleSurfaceForce(new double[] {0.0, 0.0, -10000}, BVPType.COMMON, RHS, pressBdLabel);
		fs.applyBC_MBN(A, Direct.All, fixBdLabel);
		fs.applyBC_MBN(RHS, Direct.All, 0.0, fixBdLabel);
		NewIterSSolver solver = new NewIterSSolver(A);
		solver.PCGSSOR(RHS, x, 1.65, 2);
		
		double[][] disp = new double[3][mesh3D.getNv()];
		
		fs.extractComponent(x, disp);
		mesh3D.toTecplot("matlabResults\\BracketDeformation.dat", disp);

		//mesh3.toTecplot("Elasticity3D.dat", uvw);
		System.out.println("u min: " + MVO.min(disp[0]) + "\tmax: " + MVO.max(disp[0]));
		System.out.println("v min: " + MVO.min(disp[1]) + "\tmax: " + MVO.max(disp[1]));
		System.out.println("w min: " + MVO.min(disp[2]) + "\tmax: " + MVO.max(disp[2]));
		*/
		
	}

}
